
> ## Load objects
> stayer_coef <- rio::import("est/stayer_analyses_model_coefficients_main.csv") %>%
+   dplyr::filter(model_id %in% c(209, 224)) %>%
+   dplyr::select(-V1)

> stayer_desc <- rio::import("est/stayer_analyses_descriptives_main.csv") %>%
+   dplyr::filter(model_id %in% c(209, 224))

> stayer_info <- rio::import("est/stayer_analyses_model_info_main.csv") %>%
+   dplyr::filter(model_id %in% c(209, 224))

> ## Harmonize variable names
> colnames(stayer_coef) <- colnames(stayer_coef) %>%
+   gsub("X.", "", .) %>%
+   gsub("\\._", "\\_", .) %>%
+   gsub(".cold_rent_sqm_cwu_imp1.1", ".1.cold_rent_sqm_cwu_imp1", .) %>%
+   gsub(".cold_rent_sqm_cwu_imp2.1", ".1.cold_rent_sqm_cwu_imp2", .) %>%
+   gsub(".cold_rent_sqm_cwu_imp3.1", ".1.cold_rent_sqm_cwu_imp3", .) %>%
+   gsub(".cold_rent_sqm_cwu_imp4.1", ".1.cold_rent_sqm_cwu_imp4", .) %>%
+   gsub(".cold_rent_sqm_cwu_imp5.1", ".1.cold_rent_sqm_cwu_imp5", .) %>%
+   gsub(".cold_rent_sqm_umn_imp1.1", ".1.cold_rent_sqm_umn_imp1", .) %>%
+   gsub(".cold_rent_sqm_umn_imp2.1", ".1.cold_rent_sqm_umn_imp2", .) %>%
+   gsub(".cold_rent_sqm_umn_imp3.1", ".1.cold_rent_sqm_umn_imp3", .) %>%
+   gsub(".cold_rent_sqm_umn_imp4.1", ".1.cold_rent_sqm_umn_imp4", .) %>%
+   gsub(".cold_rent_sqm_umn_imp5.1", ".1.cold_rent_sqm_umn_imp5", .) %>%
+   gsub("gtyp5urban...500k_imp1.1", "gtyp5urban...500k.1_imp1", .) %>%
+   gsub("gtyp5urban...500k_imp2.1", "gtyp5urban...500k.1_imp2", .) %>%
+   gsub("gtyp5urban...500k_imp3.1", "gtyp5urban...500k.1_imp3", .) %>%
+   gsub("gtyp5urban...500k_imp4.1", "gtyp5urban...500k.1_imp4", .) %>%
+   gsub("gtyp5urban...500k_imp5.1", "gtyp5urban...500k.1_imp5", .)

> ## Recode x_names in stayer_coef (to match column names)
> stayer_coef <- stayer_coef %>%
+   mutate(x_names = x_names %>%
+            as.character() %>%
+            ifelse(
+              startsWith(., "gtyp5urban > 500k"),
+              paste0(substr(., 1, nchar("gtyp5urban > 500k")),
+                     ".1",
+                     substr(., nchar("gtyp5urban > 500k") + 1, nchar(.))),
+              .
+            )) %>%
+   mutate(x_names = x_names %>%
+            as.character() %>%
+            ifelse(endsWith(., "gtyp5urban > 500k"),
+                   paste0(., ".1"),
+                   .)) %>%
+   mutate(
+     x_names = x_names %>%
+       gsub("\\(", "", .) %>%
+       gsub("\\)", "", .) %>%
+       gsub("<", "\\.", .) %>%
+       gsub(">", "\\.", .) %>%
+       gsub("=", "\\.", .) %>%
+       gsub(" ", "\\.", .) %>%
+       gsub("\\:", "\\.", .) %>%
+       gsub("&", "\\.", .) %>%
+       gsub("-", "\\.", .) %>%
+       gsub("\\`", "", .)
+   )

> ## Split stayer_coef by model_id, drop empty columns
> stayer_coef <- stayer_coef %>%
+   group_by(model_id) %>%
+   group_split %>%
+   lapply(select_if, all_na)

> ## Stack by imputations
> stayer_est <- stayer_coef %>%
+   lapply(.,
+          function(mod) {
+            varying <- mod %>%
+              dplyr::select(contains("_imp")) %>%
+              names() %>%
+              split(
+                .,
+                mod %>%
+                  dplyr::select(contains("_imp")) %>%
+                  names() %>%
+                  sapply(function(name)
+                    substr(name, 1, nchar(name) - 5)) %>%
+                  factor(
+                    .,
+                    levels = mod %>%
+                      dplyr::select(contains("_imp")) %>%
+                      names() %>%
+                      sapply(function(name)
+                        substr(name, 1, nchar(name) - 5)) %>%
+                      unique()
+                  )
+              ) %>%
+              lapply(function(names)
+                sapply(names, function(name)
+                  which(names(mod) == name)))
+            mod <- mod %>%
+              as.data.frame() %>%
+              reshape(
+                direction = "long",
+                varying = varying,
+                v.names = names(varying),
+                timevar = "imp",
+                times = 1:5
+              ) %>%
+              dplyr::select(-id)
+            
+            return(mod)
+          })

> ## Split by model_id and imp
> stayer_est <- stayer_est %>%
+   lapply(.,
+          function(mod) {
+            mod %>%
+              group_by(imp) %>%
+              group_split()
+          })

> ## Simulate and combine
> stayer_imp_sim <- stayer_est %>%
+   lapply(.,
+          function (mod) {
+            b_sim <- lapply(mod,
+                            function (imp) {
+                              x_names <- imp$x_names
+                              b <- as.matrix(imp[, 4])
+                              vcov <- as.matrix(imp[, 5:ncol(imp)])
+                              vcov <- vcov[, x_names]
+                              set.seed(20230329)
+                              sim <- MASS::mvrnorm(10000L, b, vcov)
+                              
+                              return(list(x_names = x_names,
+                                          sim = sim))
+                            }) 
+            x_names <- b_sim[[1]]$x_names
+            b_sim <- lapply(b_sim, function (obj)
+              obj$sim) %>%
+              do.call(rbind, .) %>%
+              apply(., 2, quantile, c(.5, .025, .975)) %>%
+              round(3) %>%
+              format(nbig = 2, nsmall = 3) %>%
+              apply(., 2, function (x) {
+                paste0("\\makecell{", "$", x[1], "$ \\\\ $", "[", x[2], ",", x[3], "]$", "}")
+              }) %>% as.matrix()
+            b_sim <- cbind.data.frame(x_names,
+                                      b_sim)
+          })

> dict <- c(
+   "cmr_arm_cwu" = "Local market rent (EUR/sqm)",
+   "log_hinc_eq_cwu" = "Equiv. household income (log)",
+   "prop_personal_hinc_cwu" = "Proportion personal income",
+   "hh_prop_ecact_cwu" = "Proportion econ. active household members",
+   "hh_mmb_cwu" = "Number of household members",
+   "lm_part_Atypical_cwu" = "\\hspace{5mm}   Atypical employment",
+   "lm_part_Inactive_cwu" = "\\hspace{5mm}   Economically inactive",
+   "lm_part_Unemployed_cwu" = "\\hspace{5mm}   Unemployed",
+   "lm_part_In_Education_cwu" = "\\hspace{5mm}   In education",
+   "lm_part_Pensioners_cwu" = "\\hspace{5mm}   Retired",
+   "cold_rent_sqm_cwu" = "Household rent (EUR/sqm)",
+   "cmr_arm_cwu.log_hinc_eq" = "Local market rent (demeaned) $\\times$ household income",
+   "log_hinc_eq.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income",
+   "Intercept" = "Intercept",
+   "log_hinc_eq" = "Equiv. household income (log)",
+   "age" = "Age",
+   "east" = "\\hspace{5mm}   East",
+   "fem" = "\\hspace{5mm}   Female",
+   "edu5Upper.Secondary" = "\\hspace{5mm}   Upper-secondary",
+   "edu5Post.Secondary.Non.Tertiary" = "\\hspace{5mm}   Post-secondary non-tertiary",
+   "edu5Higher.Vocational" = "\\hspace{5mm}   Higher vocational",
+   "edu5Higher.Education" = "\\hspace{5mm}   Tertiary",
+   "gtyp3suburban" = "\\hspace{5mm}   Suburban",
+   "gtyp3urban" = "\\hspace{5mm}   Urban",
+   "year_fac2015" = "\\hspace{5mm}   2015",
+   "year_fac2016" = "\\hspace{5mm}   2016",
+   "year_fac2017" = "\\hspace{5mm}   2017",
+   "year_fac2018" = "\\hspace{5mm}   2018",
+   "cmr_arm_umn" = "Local market rent (EUR/sqm)",
+   "prop_personal_hinc_umn" = "Proportion personal income",
+   "hh_prop_ecact_umn" = "Proportion econ. active household members",
+   "hh_mmb_umn" = "Number of household members",
+   "lm_part_Atypical_umn" = "\\hspace{5mm}   Atypical employment",
+   "lm_part_Inactive_umn" = "\\hspace{5mm}   Economically inactive",
+   "lm_part_Unemployed_umn" = "\\hspace{5mm}   Unemployed",
+   "lm_part_In_Education_umn" = "\\hspace{5mm}   In education",
+   "lm_part_Pensioners_umn" = "\\hspace{5mm}   Retired",
+   "cold_rent_sqm_umn" = "Household rent (EUR/sqm)",
+   "log_hinc_eq.cmr_arm_umn" = "Local market rent (mean) $\\times$ household income",
+   "log_hinc_eq.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income"
+ )

> stayer_imp_sim <- stayer_imp_sim %>%
+   lapply(
+     .,
+     . %>%
+       dplyr::mutate(
+         type = dplyr::case_when(
+           grepl("_cwu", x_names) ~ "within",
+           grepl("_umn", x_names) | x_names == "log_hinc_eq" ~ "between",
+           TRUE ~ "mixed"
+         ),
+         type = factor(
+           type,
+           levels = c("within", "mixed", "between"),
+           labels = c("within", "mixed", "between"),
+           ordered = TRUE
+         )
+       ) %>%
+       dplyr::arrange(type)
+   )

> ## Combined table (upper part)
> upper_table <- stayer_imp_sim[[1]]

> upper_table <- cbind(
+   upper_table,
+   NA_character_
+ )

> which_rows <- match(stayer_imp_sim[[2]][, "x_names"], upper_table[, "x_names"])

> upper_table[which_rows, ncol(upper_table)] <- stayer_imp_sim[[2]][, "b_sim"]

> ## Separate by effect type
> make_dummy_table <- function(table, dummy_pos, dummy_labels) {
+   if (dummy_pos[1] == 1) {
+     new_table <- rbind(c(
+       dummy_labels[1],
+       rep("", ncol(upper_table) - 1)
+     ),
+     table[dummy_pos[1]:(dummy_pos[2] - 1),])
+   } else {
+     new_table <- rbind(
+       table[1:(dummy_pos[1] - 1),],
+       c(
+         dummy_labels[1],
+         rep("", ncol(upper_table) - 1)
+       ),
+       table[dummy_pos[1]:(dummy_pos[2] - 1),]
+     )
+   }
+   dummy_pos <- c(dummy_pos, nrow(table) + 1)
+   for (n in 2:length(dummy_labels)) {
+     new_table <- rbind(
+       new_table,
+       c(
+         dummy_labels[n],
+         rep("", ncol(upper_table) - 1)
+       ),
+       table[dummy_pos[n]:(dummy_pos[n + 1] - 1), ]
+     )
+   }
+   return(new_table)
+ }

> row_pos <- c(
+   which(upper_table[, "type"] == "within")[1],
+   which(upper_table[, "type"] == "mixed")[1],
+   which(upper_table[, "type"] == "between")[1]
+ )

> row_labels <- c(
+   "\\emph{Within effects}",
+   "\\emph{Mixed effects}",
+   "\\emph{Between effects}"
+ )

> upper_table <- make_dummy_table(
+   upper_table,
+   row_pos,
+   row_labels
+ )

> upper_table <- upper_table[, c("x_names", "b_sim", "NA_character_")]

> colnames(upper_table) <- c("Variables", "Renter model", "Owner model")

> ## Dummy sets
> dummy_pos <-
+   c(
+     which(upper_table[, "Variables"] == "lm_part_Atypical_cwu"),
+     which(upper_table[, "Variables"] == "east"),
+     which(upper_table[, "Variables"] == "fem"),
+     which(upper_table[, "Variables"] == "edu5Upper.Secondary"),
+     which(upper_table[, "Variables"] == "gtyp3suburban"),
+     which(upper_table[, "Variables"] == "year_fac2015"),
+     which(upper_table[, "Variables"] == "lm_part_Atypical_umn")
+   )

> dummy_labels <- c(
+   "\\emph{Labor market status (ref: Full-time employment)}",
+   "\\emph{East/West residence (ref: West)}",
+   "\\emph{Sex (ref: Male)}",
+   "\\emph{Education (ref: $<$ Upper secondary)}",
+   "\\emph{Locality (ref: Rural)}",
+   "\\emph{Year (ref: 2014)}",
+   "\\emph{Labor market status (ref: Full-time employment)}"
+ )

> ## Map variable names
> upper_table[, "Variables"] <-
+   sapply(upper_table[, "Variables"], function(name)
+     ifelse(name %in% names(dict),
+            dict[which(names(dict) == name)],
+            name))

> upper_table <- make_dummy_table(
+   upper_table,
+   dummy_pos,
+   dummy_labels
+ )

> ## Combined table (lower part, num obs)
> lower_table_n <- stayer_info %>%
+   dplyr::select(starts_with("N_")) %>%
+   dplyr::mutate_all(as.character) %>%
+   t()

> lower_table_n <- cbind(
+   c(
+     "Observations",
+     "Individuals",
+     "Postcode areas",
+     "Counties/cities"
+   ),
+   lower_table_n
+ )

> ## Combined table (lower part, variance terms)
> lower_table_s <- stayer_info %>%
+   dplyr::select(model_id, starts_with("sigma_")) %>%
+   tidyr::pivot_longer(-model_id) %>%
+   dplyr::mutate(imp = dplyr::case_when(
+     grepl(".*_imp1", name) ~ 1L,
+     grepl(".*_imp2", name) ~ 2L,
+     grepl(".*_imp3", name) ~ 3L,
+     grepl(".*_imp4", name) ~ 4L,
+     grepl(".*_imp5", name) ~ 5L)
+   ) %>%
+   dplyr::mutate(name = substr(name, 1, nchar(name) - 5L)) %>%
+   dplyr::group_by (model_id, name) %>%
+   dplyr::summarize(value = mean(value)) %>%
+   dplyr::mutate(
+     name = dplyr::case_when(
+       name == "sigma_res" ~ "$\\sigma_{\\text{Observations}}$",
+       name == "sigma_id" ~ "$\\sigma_{\\text{Individuals}}$",
+       name == "sigma_plz" ~ "$\\sigma_{\\text{Postcode areas}}$",
+       name == "sigma_kkz" ~ "$\\sigma_{\\text{Counties/cities}}$"
+     ),
+     value = format(round(value, 3), nbig = 2, nsmall = 3)
+   ) %>%
+   tidyr::pivot_wider(names_from = model_id) %>%
+   as.matrix()

> ## Combine
> combined_table <- rbind(
+   as.matrix(upper_table),
+   c(
+     "$N$",
+     "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
+     rep("", ncol(upper_table) - 2)
+   ),
+   lower_table_n,
+   c(
+     "\\emph{Standard deviations of intercepts}",
+     "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
+     rep("", ncol(upper_table) - 2)
+   ),
+   lower_table_s
+ )

> ## TeX output
> print(
+   xtable(
+     combined_table,
+     align = c("l", "l", "c", "c"),
+     label = "tab:reg1",
+     caption = paste0(
+       "Coefficients and simulation-based 95\\% confidence intervals from hierarchical linear within-between models, estimated across $M=5$ imputations."
+     )
+   ),
+   include.rownames = FALSE,
+   include.colnames = TRUE,
+   sanitize.text.function = identity,
+   type = "latex",
+   hline.after = c(-1,
+                   -1,
+                   0,
+                   which(combined_table[, "Variables"] == "\\emph{Mixed effects}") - 1,
+                   which(combined_table[, "Variables"] == "\\emph{Between effects}") - 1,
+                   which(combined_table[, "Variables"] == "$N$") - 1,
+                   nrow(combined_table),
+                   nrow(combined_table)),
+   size = "\\small",
+   comment = TRUE,
+   floating = FALSE,
+   tabular.environment = "longtable",
+   file = "tex/reg1_2024.tex"
+ )

> ## Load objects
> stayer_coef <- rio::import("est/stayer_analyses_model_coefficients_main.csv") %>%
+   dplyr::filter(model_id %in% c(254, 269)) %>%
+   dplyr::select(-V1)

> stayer_desc <- rio::import("est/stayer_analyses_descriptives_main.csv") %>%
+   dplyr::filter(model_id %in% c(254, 269))

> stayer_info <- rio::import("est/stayer_analyses_model_info_main.csv") %>%
+   dplyr::filter(model_id %in% c(254, 269))

> ## Harmonize variable names
> colnames(stayer_coef) <- colnames(stayer_coef) %>%
+   gsub("X.", "", .) %>%
+   gsub("\\._", "\\_", .) %>%
+   gsub(".cold_rent_sqm_cwu_imp1.1", ".1.cold_rent_sqm_cwu_imp1", .) %>%
+   gsub(".cold_rent_sqm_cwu_imp2.1", ".1.cold_rent_sqm_cwu_imp2", .) %>%
+   gsub(".cold_rent_sqm_cwu_imp3.1", ".1.cold_rent_sqm_cwu_imp3", .) %>%
+   gsub(".cold_rent_sqm_cwu_imp4.1", ".1.cold_rent_sqm_cwu_imp4", .) %>%
+   gsub(".cold_rent_sqm_cwu_imp5.1", ".1.cold_rent_sqm_cwu_imp5", .) %>%
+   gsub(".cold_rent_sqm_umn_imp1.1", ".1.cold_rent_sqm_umn_imp1", .) %>%
+   gsub(".cold_rent_sqm_umn_imp2.1", ".1.cold_rent_sqm_umn_imp2", .) %>%
+   gsub(".cold_rent_sqm_umn_imp3.1", ".1.cold_rent_sqm_umn_imp3", .) %>%
+   gsub(".cold_rent_sqm_umn_imp4.1", ".1.cold_rent_sqm_umn_imp4", .) %>%
+   gsub(".cold_rent_sqm_umn_imp5.1", ".1.cold_rent_sqm_umn_imp5", .) %>%
+   gsub("gtyp5urban...500k_imp1.1", "gtyp5urban...500k.1_imp1", .) %>%
+   gsub("gtyp5urban...500k_imp2.1", "gtyp5urban...500k.1_imp2", .) %>%
+   gsub("gtyp5urban...500k_imp3.1", "gtyp5urban...500k.1_imp3", .) %>%
+   gsub("gtyp5urban...500k_imp4.1", "gtyp5urban...500k.1_imp4", .) %>%
+   gsub("gtyp5urban...500k_imp5.1", "gtyp5urban...500k.1_imp5", .)

> ## Recode x_names in stayer_coef (to match column names)
> stayer_coef <- stayer_coef %>%
+   mutate(x_names = x_names %>%
+            as.character() %>%
+            ifelse(
+              startsWith(., "gtyp5urban > 500k"),
+              paste0(substr(., 1, nchar("gtyp5urban > 500k")),
+                     ".1",
+                     substr(., nchar("gtyp5urban > 500k") + 1, nchar(.))),
+              .
+            )) %>%
+   mutate(x_names = x_names %>%
+            as.character() %>%
+            ifelse(endsWith(., "gtyp5urban > 500k"),
+                   paste0(., ".1"),
+                   .)) %>%
+   mutate(
+     x_names = x_names %>%
+       gsub("\\(", "", .) %>%
+       gsub("\\)", "", .) %>%
+       gsub("<", "\\.", .) %>%
+       gsub(">", "\\.", .) %>%
+       gsub("=", "\\.", .) %>%
+       gsub(" ", "\\.", .) %>%
+       gsub("\\:", "\\.", .) %>%
+       gsub("&", "\\.", .) %>%
+       gsub("-", "\\.", .) %>%
+       gsub("\\`", "", .)
+   )

> ## Split stayer_coef by model_id, drop empty columns
> stayer_coef <- stayer_coef %>%
+   group_by(model_id) %>%
+   group_split %>%
+   lapply(select_if, all_na)

> ## Stack by imputations
> stayer_est <- stayer_coef %>%
+   lapply(.,
+          function(mod) {
+            varying <- mod %>%
+              dplyr::select(contains("_imp")) %>%
+              names() %>%
+              split(
+                .,
+                mod %>%
+                  dplyr::select(contains("_imp")) %>%
+                  names() %>%
+                  sapply(function(name)
+                    substr(name, 1, nchar(name) - 5)) %>%
+                  factor(
+                    .,
+                    levels = mod %>%
+                      dplyr::select(contains("_imp")) %>%
+                      names() %>%
+                      sapply(function(name)
+                        substr(name, 1, nchar(name) - 5)) %>%
+                      unique()
+                  )
+              ) %>%
+              lapply(function(names)
+                sapply(names, function(name)
+                  which(names(mod) == name)))
+            mod <- mod %>%
+              as.data.frame() %>%
+              reshape(
+                direction = "long",
+                varying = varying,
+                v.names = names(varying),
+                timevar = "imp",
+                times = 1:5
+              ) %>%
+              dplyr::select(-id)
+            
+            return(mod)
+          })

> ## Split by model_id and imp
> stayer_est <- stayer_est %>%
+   lapply(.,
+          function(mod) {
+            mod %>%
+              group_by(imp) %>%
+              group_split()
+          })

> ## Simulate and combine
> stayer_imp_sim <- stayer_est %>%
+   lapply(.,
+          function (mod) {
+            b_sim <- lapply(mod,
+                            function (imp) {
+                              x_names <- imp$x_names
+                              b <- as.matrix(imp[, 4])
+                              vcov <- as.matrix(imp[, 5:ncol(imp)])
+                              vcov <- vcov[, x_names]
+                              set.seed(20230329)
+                              sim <- MASS::mvrnorm(10000L, b, vcov)
+                              
+                              return(list(x_names = x_names,
+                                          sim = sim))
+                            }) 
+            x_names <- b_sim[[1]]$x_names
+            b_sim <- lapply(b_sim, function (obj)
+              obj$sim) %>%
+              do.call(rbind, .) %>%
+              apply(., 2, quantile, c(.5, .025, .975)) %>%
+              round(3) %>%
+              format(nbig = 2, nsmall = 3) %>%
+              apply(., 2, function (x) {
+                paste0("\\makecell{", "$", x[1], "$ \\\\ $", "[", x[2], ",", x[3], "]$", "}")
+              }) %>% as.matrix()
+            b_sim <- cbind.data.frame(x_names,
+                                      b_sim)
+          })

> dict <- c(
+   "cmr_arm_cwu" = "Local market rent (EUR/sqm)",
+   "log_hinc_eq_cwu" = "Equiv. household income (log)",
+   "prop_personal_hinc_cwu" = "Proportion personal income",
+   "hh_prop_ecact_cwu" = "Proportion econ. active household members",
+   "hh_mmb_cwu" = "Number of household members",
+   "lm_part_Atypical_cwu" = "\\hspace{5mm}   Atypical employment",
+   "lm_part_Inactive_cwu" = "\\hspace{5mm}   Economically inactive",
+   "lm_part_Unemployed_cwu" = "\\hspace{5mm}   Unemployed",
+   "lm_part_In_Education_cwu" = "\\hspace{5mm}   In education",
+   "lm_part_Pensioners_cwu" = "\\hspace{5mm}   Retired",
+   "cold_rent_sqm_cwu" = "Household rent (EUR/sqm)",
+   "cmr_arm_cwu.log_hinc_eq" = "Local market rent (demeaned) $\\times$ household income",
+   "log_hinc_eq.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income",
+   "Intercept" = "Intercept",
+   "log_hinc_eq" = "Equiv. household income (log)",
+   "age" = "Age",
+   "east" = "\\hspace{5mm}   East",
+   "fem" = "\\hspace{5mm}   Female",
+   "edu5Upper.Secondary" = "\\hspace{5mm}   Upper-secondary",
+   "edu5Post.Secondary.Non.Tertiary" = "\\hspace{5mm}   Post-secondary non-tertiary",
+   "edu5Higher.Vocational" = "\\hspace{5mm}   Higher vocational",
+   "edu5Higher.Education" = "\\hspace{5mm}   Tertiary",
+   "gtyp3suburban" = "\\hspace{5mm}   Suburban",
+   "gtyp3urban" = "\\hspace{5mm}   Urban",
+   "year_fac2015" = "\\hspace{5mm}   2015",
+   "year_fac2016" = "\\hspace{5mm}   2016",
+   "year_fac2017" = "\\hspace{5mm}   2017",
+   "year_fac2018" = "\\hspace{5mm}   2018",
+   "cmr_arm_umn" = "Local market rent (EUR/sqm)",
+   "prop_personal_hinc_umn" = "Proportion personal income",
+   "hh_prop_ecact_umn" = "Proportion econ. active household members",
+   "hh_mmb_umn" = "Number of household members",
+   "lm_part_Atypical_umn" = "\\hspace{5mm}   Atypical employment",
+   "lm_part_Inactive_umn" = "\\hspace{5mm}   Economically inactive",
+   "lm_part_Unemployed_umn" = "\\hspace{5mm}   Unemployed",
+   "lm_part_In_Education_umn" = "\\hspace{5mm}   In education",
+   "lm_part_Pensioners_umn" = "\\hspace{5mm}   Retired",
+   "cold_rent_sqm_umn" = "Household rent (EUR/sqm)",
+   "log_hinc_eq.cmr_arm_umn" = "Local market rent (mean) $\\times$ household income",
+   "log_hinc_eq.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income",
+   "cmr_arm_cwu.gtyp3suburban" = "Market rent (demeaned)  $\\times$ suburban locality",
+   "cmr_arm_cwu.gtyp3urban" = "Market rent (demeaned) $\\times$ urban locality",
+   "log_hinc_eq.gtyp3suburban" = "Household income $\\times$ suburban locality",
+   "log_hinc_eq.gtyp3urban" = "Household income $\\times$ urban locality",
+   "gtyp3suburban.cmr_arm_umn" = "Market rent (mean)  $\\times$ suburban locality",
+   "gtyp3urban.cmr_arm_umn" = "Market rent (mean) $\\times$ urban locality",
+   "gtyp3suburban.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income $\\times$ suburban locality",
+   "gtyp3urban.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income $\\times$ urban locality",
+   "gtyp3suburban.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income $\\times$ suburban locality",
+   "gtyp3urban.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income $\\times$ urban locality",
+   "cmr_arm_cwu.log_hinc_eq.gtyp3suburban" = "Market rent (demeaned) $\\times$ household income $\\times$ suburban locality",
+   "cmr_arm_cwu.log_hinc_eq.gtyp3urban" = "Market rent (demeaned) $\\times$ household income $\\times$ urban locality",
+   "log_hinc_eq.gtyp3suburban.cmr_arm_umn" = "Market rent (mean) $\\times$ household income $\\times$ suburban locality",
+   "log_hinc_eq.gtyp3urban.cmr_arm_umn" = "Market rent (mean) $\\times$ household income $\\times$ urban locality",
+   "log_hinc_eq.gtyp3suburban.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income $\\times$ suburban locality",
+   "log_hinc_eq.gtyp3urban.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income $\\times$ urban locality",
+   "log_hinc_eq.gtyp3suburban.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income $\\times$ suburban locality",
+   "log_hinc_eq.gtyp3urban.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income $\\times$ urban locality" 
+ )

> stayer_imp_sim <- stayer_imp_sim %>%
+   lapply(
+     .,
+     . %>%
+       dplyr::mutate(
+         type = dplyr::case_when(
+           grepl("_cwu", x_names) ~ "within",
+           grepl("_umn", x_names) | x_names == "log_hinc_eq" ~ "between",
+           TRUE ~ "mixed"
+         ),
+         type = factor(
+           type,
+           levels = c("within", "mixed", "between"),
+           labels = c("within", "mixed", "between"),
+           ordered = TRUE
+         )
+       ) %>%
+       dplyr::arrange(type)
+   )

> ## Combined table (upper part)
> upper_table2 <- stayer_imp_sim[[1]]

> upper_table2 <- cbind(
+   upper_table2,
+   NA_character_
+ )

> which_rows <- match(stayer_imp_sim[[2]][, "x_names"], upper_table2[, "x_names"])

> upper_table2[which_rows, ncol(upper_table2)] <- stayer_imp_sim[[2]][, "b_sim"]

> ## Separate by effect type
> make_dummy_table <- function(table, dummy_pos, dummy_labels) {
+   if (dummy_pos[1] == 1) {
+     new_table <- rbind(c(
+       dummy_labels[1],
+       rep("", ncol(upper_table) - 1)
+     ),
+     table[dummy_pos[1]:(dummy_pos[2] - 1),])
+   } else {
+     new_table <- rbind(
+       table[1:(dummy_pos[1] - 1),],
+       c(
+         dummy_labels[1],
+         rep("", ncol(upper_table) - 1)
+       ),
+       table[dummy_pos[1]:(dummy_pos[2] - 1),]
+     )
+   }
+   dummy_pos <- c(dummy_pos, nrow(table) + 1)
+   for (n in 2:length(dummy_labels)) {
+     new_table <- rbind(
+       new_table,
+       c(
+         dummy_labels[n],
+         rep("", ncol(upper_table) - 1)
+       ),
+       table[dummy_pos[n]:(dummy_pos[n + 1] - 1), ]
+     )
+   }
+   return(new_table)
+ }

> row_pos <- c(
+   which(upper_table2[, "type"] == "within")[1],
+   which(upper_table2[, "type"] == "mixed")[1],
+   which(upper_table2[, "type"] == "between")[1]
+ )

> row_labels <- c(
+   "\\emph{Within effects}",
+   "\\emph{Mixed effects}",
+   "\\emph{Between effects}"
+ )

> upper_table2 <- make_dummy_table(
+   upper_table2,
+   row_pos,
+   row_labels
+ )

> upper_table2 <- upper_table2[, c("x_names", "b_sim", "NA_character_")]

> colnames(upper_table2) <- c("Variables", "Renter model", "Owner model")

> ## Dummy sets
> dummy_pos <-
+   c(
+     which(upper_table2[, "Variables"] == "lm_part_Atypical_cwu"),
+     which(upper_table2[, "Variables"] == "east"),
+     which(upper_table2[, "Variables"] == "fem"),
+     which(upper_table2[, "Variables"] == "edu5Upper.Secondary"),
+     which(upper_table2[, "Variables"] == "gtyp3suburban"),
+     which(upper_table2[, "Variables"] == "year_fac2015"),
+     which(upper_table2[, "Variables"] == "lm_part_Atypical_umn")
+   )

> dummy_labels <- c(
+   "\\emph{Labor market status (ref: Full-time employment)}",
+   "\\emph{East/West residence (ref: West)}",
+   "\\emph{Sex (ref: Male)}",
+   "\\emph{Education (ref: $<$ Upper secondary)}",
+   "\\emph{Locality (ref: Rural)}",
+   "\\emph{Year (ref: 2014)}",
+   "\\emph{Labor market status (ref: Full-time employment)}"
+ )

> dummy_labels <- dummy_labels[order(dummy_pos)]

> dummy_pos <- dummy_pos[order(dummy_pos)]

> ## Map variable names
> upper_table2[, "Variables"] <-
+   sapply(upper_table2[, "Variables"], function(name)
+     ifelse(name %in% names(dict),
+            dict[which(names(dict) == name)],
+            name))

> upper_table2 <- make_dummy_table(
+   upper_table2,
+   dummy_pos,
+   dummy_labels
+ )

> ## Combined table (lower part, num obs)
> lower_table2_n <- stayer_info %>%
+   dplyr::select(starts_with("N_")) %>%
+   dplyr::mutate_all(as.character) %>%
+   t()

> lower_table2_n <- cbind(
+   c(
+     "Observations",
+     "Individuals",
+     "Postcode areas",
+     "Counties/cities"
+   ),
+   lower_table2_n
+ )

> ## Combined table (lower part, variance terms)
> lower_table2_s <- stayer_info %>%
+   dplyr::select(model_id, starts_with("sigma_")) %>%
+   tidyr::pivot_longer(-model_id) %>%
+   dplyr::mutate(imp = dplyr::case_when(
+     grepl(".*_imp1", name) ~ 1L,
+     grepl(".*_imp2", name) ~ 2L,
+     grepl(".*_imp3", name) ~ 3L,
+     grepl(".*_imp4", name) ~ 4L,
+     grepl(".*_imp5", name) ~ 5L)
+   ) %>%
+   dplyr::mutate(name = substr(name, 1, nchar(name) - 5L)) %>%
+   dplyr::group_by (model_id, name) %>%
+   dplyr::summarize(value = mean(value)) %>%
+   dplyr::mutate(
+     name = dplyr::case_when(
+       name == "sigma_res" ~ "$\\sigma_{\\text{Observations}}$",
+       name == "sigma_id" ~ "$\\sigma_{\\text{Individuals}}$",
+       name == "sigma_plz" ~ "$\\sigma_{\\text{Postcode areas}}$",
+       name == "sigma_kkz" ~ "$\\sigma_{\\text{Counties/cities}}$"
+     ),
+     value = format(round(value, 3), nbig = 2, nsmall = 3)
+   ) %>%
+   tidyr::pivot_wider(names_from = model_id) %>%
+   as.matrix()

> ## Combine
> combined_table2 <- rbind(
+   as.matrix(upper_table2),
+   c(
+     "$N$",
+     "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
+     rep("", ncol(upper_table2) - 2)
+   ),
+   lower_table2_n,
+   c(
+     "\\emph{Standard deviations of intercepts}",
+     "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
+     rep("", ncol(upper_table2) - 2)
+   ),
+   lower_table2_s
+ )

> ## TeX output
> print(
+   xtable(
+     combined_table,
+     align = c("l", "l", "c", "c"),
+     label = "tab:reg2",
+     caption = paste0(
+       "Coefficients and simulation-based 95\\% confidence intervals from hierarchical linear within-between models, estimated across $M=5$ imputations."
+     )
+   ),
+   include.rownames = FALSE,
+   include.colnames = TRUE,
+   sanitize.text.function = identity,
+   type = "latex",
+   hline.after = c(-1,
+                   -1,
+                   0,
+                   which(combined_table[, "Variables"] == "\\emph{Mixed effects}") - 1,
+                   which(combined_table[, "Variables"] == "\\emph{Between effects}") - 1,
+                   which(combined_table[, "Variables"] == "$N$") - 1,
+                   nrow(combined_table),
+                   nrow(combined_table)),
+   size = "\\small",
+   comment = TRUE,
+   floating = FALSE,
+   tabular.environment = "longtable",
+   file = "tex/reg2_2024.tex"
+ )

> ## Reduced table
> rows <- c(
+   "Local market rent (EUR/sqm)",
+   "Local market rent (demeaned) $\\times$ household income",
+   "Local market rent (demeaned) $\\times$ household income",
+   "Market rent (demeaned)  $\\times$ suburban locality",
+   "Market rent (demeaned) $\\times$ urban locality",
+   "Market rent (demeaned) $\\times$ household income $\\times$ suburban locality",
+   "Market rent (demeaned) $\\times$ household income $\\times$ urban locality"
+ )

> upper_table3 <- upper_table2 %>%
+   dplyr::slice(1:which(Variables == "\\emph{Mixed effects}")) %>%
+   dplyr::filter(Variables %in% rows) %>%
+   dplyr::left_join(upper_table %>%
+                      dplyr::slice(1:which(Variables == "\\emph{Mixed effects}")) %>%
+                      dplyr::filter(Variables %in% rows),
+                    by = "Variables") %>%
+   dplyr::mutate(
+     Variables = ifelse(
+       Variables == "Local market rent (EUR/sqm)",
+       "Local market rent (demeaned)",
+       Variables
+     )
+   ) %>%
+   dplyr::select(
+     Variables,
+     `Renter model.y`,
+     `Owner model.y`,
+     `Renter model.x`,
+     `Owner model.x`
+   ) %>%
+   dplyr::rename(
+     `\\makecell{Renter model \\\\ (Overall)}` = `Renter model.y`,
+     `\\makecell{ Owner model \\\\ (Overall)}` = `Owner model.y`,
+     `\\makecell{Renter model \\\\ (Locality-specific)}` = `Renter model.x`,
+     `\\makecell{ Owner model \\\\ (Locality-specific)}` = `Owner model.x`
+   )

> combined_table3 <- rbind(
+   as.matrix(upper_table3),
+   c(
+     "$N$",
+     "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
+     rep("", ncol(upper_table3) - 2)
+   ),
+   cbind(lower_table_n, lower_table2_n[, 2:3]),
+   c(
+     "\\emph{Standard deviations of intercepts}",
+     "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
+     rep("", ncol(upper_table3) - 2)
+   ),
+   cbind(lower_table_s, lower_table2_s[, 2:3])
+ )

> ## TeX output
> print(
+   xtable(
+     combined_table3,
+     align = c("l", "l", "c", "c", "c", "c"),
+     label = "tab:reg3",
+     caption = paste0(
+       "Short regression table, selectively reporting constitutive terms of the marginal effects presented in Figs. \\ref{fig:afd_inc} and \\ref{fig:afd_region}. Coefficients and simulation-based 95\\% confidence intervals from hierarchical linear within-between models, estimated across $M=5$ imputations. For full regression output, see Tables \\ref{tab:reg1} and \\ref{tab:reg2} in the Online Appendix."
+     )
+   ),
+   include.rownames = FALSE,
+   include.colnames = TRUE,
+   sanitize.text.function = identity,
+   type = "latex",
+   hline.after = c(-1,
+                   -1,
+                   0,
+                   which(combined_table3[, "Variables"] == "$N$") - 1,
+                   nrow(combined_table3),
+                   nrow(combined_table3)),
+   size = "\\tiny",
+   comment = TRUE,
+   floating = FALSE,
+   tabular.environment = "longtable",
+   file = "tex/reg3_2024.tex"
+ )
